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Introduction.    The  Institute  of  Mathematical  Sciences 
at  New  York  University,  of  which  the  AEC  Computing  Facility 
Is  a  part,  has  been  authorized  by  the  AEC  to  use  some  of 
the  computing  time  set  aside  for  research.   The  following 
sections  provide  brief  descriptions  of  the  problems  which 
have  been  studied  as  part  of  this  program  during  the  past 
fiscal  year.   Since  this  Is  a  summary  report,  the  statements 
of  the  problems  are  not  complete  and  do  not  contain  much 
detail.   In  many  cases  the  computational  studies  are 
parts  of  broader  research  programs. 

Inquiries  about  any  of  the  problems  are  welcomed. 
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Problem  No.  1.12:   The  passage  of  a  shock-wave  over  an 

obstruction  In  a  channel 


Investigators: 


Peter  D.  Lax  and  Glenn  Lewis 


When  a  shock  wave  travelling  in  a  channel  encounters 
a  narrowing  or  an  obstacle  in  the  channel  (see  Figure  l), 
part  of  it  is  reflected,  part  transmitted.   What  is  of 
interest  is  to  calculate  the  strength  of  the  transmitted 
and  reflected  waves,  and  the  shape  of  the  pressure  pulse 
behind  these  waves. 


Shock-wave 
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Figure   1. 

The  calculation  of  the  effect  of  a  blast  wave  on 
structures  presents  a  similar  problem.   There  the  under- 
lying domain  is  a  half space  as  in  Figure  2: 
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Figure  2. 
The  quantity  of  Interest  here  is  the  distribution  of  pressure 
over  the  structure  as  a  function  of  time. 
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The  method  of  calculation;   Eulerian  coordinates  were 
used,  and  the  equation  of  state  was  taken  to  be  a  -y-law, 
y   =  1.4.   The  equations  of  motion  were  written  in  conserva- 
tion form: 

Pt  +  '"x  -^  "y  =  ° 

ra^  +  (um+p)^  +  ^^™V  "  ° 

"t  "^  ^^"^x  "^  (^"+P)y  =  ° 

E^  +  (u(E+p))^  +  (v(E+p))y  =  0. 

Here  m  and  n  abbreviate  momenta  In  the  x  and  y 
direction,  and  E  is  total  energy  per  unit  volume 

m  =  pu,  n  =  pv,  E  =  pe  +  2  pu  ,  e  =  ^P_j^^   . 

p,  m,  n  and  E  were  the  unknown  functions  used  In  the  compu- 
tations; the  others  were  expressed  as  functions  of  them.   A 
rectangular  grid  staggered  In  time  was  used,  i.e.  at  even 
time  cycles  all  quantities  were  expressed  at  even  lattice 
points,  at  odd  times  at  half  points  (see  Figure  3). 
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Figure  3. 

•  position  of  quantities  at  even  time  step 
X  position  of  quantities  at  odd  time  step 

Time  derivatives  were  replaced  by  averaged  forward  difference 
quotients,  i.e.  for  n  odd: 
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A        ^     n+1        Pk- 1/2 , 1-  \/2'^\^^y2l- 1/2'^  Pk- 1/2 ,  i-H/2'*"Pk+l/2 ,  i-H/g 

\P  ^  K,i — 5 — 

and  similarly  for  n  even.   Space  derivatives  were  replaced 
by  centered  difference  quotients: 


n  n  n  n 

Pv   Pk+1/2 ,  i-H/2" Pk-1/2 ,  i+1/2  .  Pk+l/2,  i-l/2"Pk-l/2 ,  l-\/2 
^  2Ax  "*"  2Ax 


At  boundary  points  reflected  values  were  used  at  the  lattice 
points  which  fell  outside  the  flow  region,  i.e.,  the  value 
of  p,  E  and  p  was  reflected  while  the  value  of  u,  v  was 
changed  so  as  to  preserve  the  tangential  flow  velocity  and 
reverse  the  normal  velocity. 

Initial  configuration:  Two  regions  with  constant 
states  connected  through  a  straight  shock,  as  shown  in 
Figure  1.   The  pressure  ratio  was  5:1. 

Stability  criterion;  A  useful  criterion  of  Friedrichs 
for  stability  was  employed,  which  for  the  scheme  under  con- 
sideration requires  that 


\\\x\    +  M.|v|    +  cy>?+|?  <   1, 


where  c  is  sound  speed,  and  \   and  p.  abbreviate 

,  _  2At        2At 

ks,   the  mesh  was  traversed  during  the  calculation,  the  current 
value  of  At  was  checked  at  each  point  and  enlarged  if 
necessary  for  the  next  time  cycle.   This  comparison  was 
arranged  so  as  to  avoid  the  extraction  of  square  roots. 

-  4  - 


Organization  of  data  in  the  machine;   The  values  of 
the  hydrodynaunlcal  variables  along  the  lattice  points  on 
one  vertical  line  (i.e.   x  =  const.)  were  stored  in  one 
block  on  tape  1.   At  any  given  time  there  were  data  from 
two  neighboring  vertical  lines  (the  left  and  the  right)  in 
the  internal  memory,  from  which  the  new  values  at  the  line  in 
between  were  computed.   These  new  values  were  dumped  on 
tape  2;  then  the  contents  of  the  right  line  were  transferred 
to  the  left  line  internally,  and  then  the  right  line  was 
filled  from  the  next  block  on  tape  1.   In  the  next  time 
cycle  the  mesh  was  swept  from  right  to  left,  reading  from 
tape  2  and  writing  on  tape  1. 

All  calculations  were  carried  out  in  fixed  decimal. 

Output;   Pressure,  density,  flow  speed  and  flow  angle 
were  tabulated  separately  to  three  significant  figures,  in 
floating  decimal  notation.   The  tables  were  replicas  of  the 
channel;  each  entry  appeared  over  the  lattice  point  whose 
value  it  represented. 

Results ;   50  time  cycles  were  ground  out,  in  about 
2  1/2  hours.   The  shock,  which  initially  started  at  some 
distance  from  the  obstacle,  was  Just  over  the  edge.   There 
was  an  Increase  in  pressure  against  the  obstacle  as  the 
reflected  shock  was  beginning  to  build  up.   The  isobars 
plotted  around  the  diffracting  corner  presented  the  kind 
of  pattern  one  would  expect  on  the  basis  of  acoustic  theory. 

The  calculation  so  far  has  given  no  inkling  of  the 
strength  of  the  reflected  and  transmitted  shocks.   We  are 
optimistic  through  that  if  run  long  enough  the  method  will 
furnish  this  information. 

-  5  - 


Problem  No.  1.13:   Underwater  explosion  bubble  oscillations  (Bubble) 
Investigators:     L.  K.  Brathwalte,  P.  Pox,  J.  B.  Keller 

When  an  explosive  detonates  underwater  it  is  immediately 
changed  into  a  gas  at  high  pressure.   This  gas,  called 
the  bubble,  expands  until  its  pressure  equals  that  of  the 
surrounding  water,  but  inertia  causes  it  to  overexpand. 
Finally  it  comes  to  rest  momentarily  at  a  lower  pressure 
but  then  the  surrounding  liquid  recompresses  it.   Then  it  expands 
again  and  this  process  repeats  a  number  of  times  before 
the  bubble  finally  breaks  up  or  vents  at  the  surface. 

The  usual  theory  of  these  oscillations  treats  the  bubble 
as  a  sphere  and  assumes  the  water  to  be  Incompressible. 
This  theory  has  been  modified  by  Dr.  J.  B.  Keller  and  Dr. 
I.  I.  Kolodner  to  take  account  of  the  compressibility  of 
the  water.   It  then  yields  damped  oscillations  of  diminish- 
ing period,  in  agreement  with  experimental  observations. 
This  modified  theory  is  described  in  "Underwater  Explosion 
Bubbles  I.   The  Effect  of  the  Compressibility  of  the  Water" 
IMM-NYU  No.  191. 

The  theory  leads  to  a  non-linear  second  order  ordinary 
differential  equation  for  the  bubble  radius  as  a  function 
of  time.   This  equation  was  programmed  and  coded  for  solu- 
tion on  Univac  by  Dr.  P.  Pox,  M.  Billings  and  L.  Brathwalte. 
A  nxomber  of  radius- time  curves  were  computed.   These  curves, 
together  with  the  theory,  will  be  published  in  the  Journal 
of  Applied  Physics  in  an  art.icle  by  J.  B.  Keller  and  I.  I. 
Kolodner  entitled  "Damping  of  Underwater  Explosion  Bubble 

Oscillations. " 
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Problem  No.  I.16:   Solution  of  elliptic  difference  equations 

(Pogo) 

Investigators:     R.  Eykaar,  K.  Gordls,  H.  B.  Keller 

A  large  class  of  finite-difference  approximations  to 
linear  elliptic  partial  differential  equations  (i.e. 
Laplace's  equation,  equations  of  neutron  diffusion,  elastic 
plate  and  membrane  equations,  etc.)  can  be  written  as  a 
system  of  linear  algebraic  equations  in  the  form 

(1)  (I  -  H)  i  =  C   . 

Here   I  is  the  unit  matrix;   H  is  a  square  matrix  whose 
elements  depend  upon  the  mesh  spacing,  the  coefficients  of 
the  original  equation,  and  the  method  of  differencing 
employed;   C  is  a  column  vector  whose  components  depend 
upon  the  mesh  spacing,  Inhomogeneous  terms  in  the  original 
equation  and  the  boundary  conditions;  and  J  is  a  column 
vector  whose  components  are  the  unknown  n\imerical  approxima- 
tion to  the  solution. 

The  order  of  this  system  is  essentially  the  number  of 
mesh  points  in  the  domain  under  consideration.   Thus  for 
problems  In  two  or  more  dimensions  the  system  (l)  may  be 
of  such  large  order  that  general  matrix  methods  of  solution 
are  impractical  or  perhaps  unacceptably  inaccurate.   However, 
the  origin  of  these  equations  leads  us  to  expect  that   H 
may  have  a  special  form  with  many  zero  elements,  and 
accordingly  some  particular  methods  for  the  solution  of 
(1)  have  been  presented.   The  Pogo  program  is  an  attempt 
T.      S.  P.  Frankel:  M.T.A.C.  4,  b^  (1950). 
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to  Investigate  and  improve  methods  for  the  nxunerical  solu- 
tion of  these  special  forms  of  the  system  (l). 

In  particular  the  class  of  methods  under  considera- 
tion are  iterative  techniques  and  may  be  generally  described 
as  follows:   The  coefficient  matrix  is  "split"  into  the 
difference  of  two  matrices 

(2)  I  -  H  =  A(a)  -  B(a)  , 

whose  elements  are  functions  of  some  real  parameter,  a,  in 

such  a  way  that  |A(a)|  /  0.   Then  from  some  Initial  guess  at  the 

solution  vector,  f^,  we  define  a  sequence  of  approximations,  ^  ,  by 

(5)  A(a)  -^   =  B(a)  J^_-^   +  C  . 

If  this  sequence  converges,  the  limit  vector  is  clearly  the 
solution  of  (l).   We  seek  "splittings"  of  the  form  (2)  and 
corresponding  values  of  the  parameter,  a,  for  which  conver- 
gence of  the  sequence  (5)  is  assured  and  is  as  rapid  as 
possible.   An  Important  restriction  to  be  Imposed  is  that 
the  resulting  equations  (5)  may  be  solved  in  some  "simple" 
straightforward  manner. 

Many  of  the  standard  iterative  procedures  for  solving 
such  systems  are  easily  cast  Into  the  above  form  and  their 
theoretical  analysis  becomes  simple  and  elegant.   These 
methods  also  suggest  quite  naturally  other  splits  which 
are  under  Investigation.   The  neighbor  relationships  in 
the  difference  equations  and  methods  of  sweeping  the  mesh 
dictated  by  efficient  use  of  digital  computers  suggest  still 
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more  splits  which  are  to  be  examined. 

For  any  split  of  form  (2)  the  theoretical  investigation 
of  convergence  requires  a  determination  of  the  eigenvalue  of 
largest  absolute  value  of  the  matrix  A"  (a)B(a).   To 
compare  rates  of  convergence  of  two  different  schemes  we 
need  only  a  comparison  of  these  eigenvalues.   However, 
for  some  splittings  which  have  been  proposed  and  used 
successfully  the  eigenvalues  are  unknown  or  have  a  compli- 
cated dependence  on  various  quantities  (i.e.  mesh  spacing, 
a,  number  of  mesh  points,  etc.)-   Such  schemes  are  being 
studied  by  means  of  numerical  experiments  in  which  con- 
vergence rates  are  observed  and  the  maximum  eigenvalues 
and  best  values  for  a  are  computed  approximately.   As  a 
basis  for  comparison  these  same  experiments  will  be  con- 
ducted for  a  scheme  (probably  extrapolated  Llebmann)  whose 
theoretical  properties  are  well  known. 

At  present  the  Pogo  program  Is  considering  second  order 
equations  in  two  dimensions  with  second  order  difference 
approximations  on  full  rectangular  meshes.   Some  recent 
investigations  by  H.  J.  Greenberg  indicate  rather  clearly 
that  much  improvement  in  efficiency  may  be  obtained  by 
considering  staggered  meshes.   The  possibilities  offered 
by  higher  order  difference  approximations  also  seem  promising. 
It  is  hoped  that  some  of  these  considerations  will  be 
included  in  future  Investigations. 
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Problem  No.  I.I8:  The  treatment  of  contact  dlscontlnutles 

In  a  Eulerlan  calculation  (Interface) 

Investigators:     R.  Buchal  and  P.  D.  Lax 


Contact  discontinuities  in  compressible  flows  behave 
like  discontinuities  of  solutions  of  linear  equations.   The 
similarity  is  present  in  finite  difference  calculation  and 
has  this  effect:  the  width  of  the  transition  region  spreads 
like  the  square  root  of  the  number  of  time  cycles,  unless 
special  cognizance  is  taken  of  the  presence  of  the  contact 
discontinuity.   If  the  contact  discontinuity  separates  two 
fluids  of  very  different  densities,  this  diffusion  across 
the  interface  causes  a  very  large  distortion.   In  Lagrange 
coordinates  it  is  possible  to  prevent  the  diffusion  of  the 
discontinuity  by  an  appropriate  differencing  of  the  equa- 
tions; this  is  not  possible  in  Euler  coordinates.   The 
reason  for  it  is  that  in  Lagrange  description  the  position 
of  the  contact  discontinuity  is  unchanged  in  the  grid;  in 
Euler  coordinates  it  is  not. 

In  our  sample  calculations  we  tried  out  some  specific 
measures  to  prevent  diffusion.   Two  steps  are  involved: 

a)  keeping  track  of  the  contact  discontinuity; 

b)  devising  special  difference  schemes  to  be 
used  at  points  whose  neighbors  straddle  the 
discontinuity  line. 

These  extra  operations  consume  negligible  time  com- 
pared to  the  main  body  of  the  calculation. 
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So  far  the  method  has  been  tried  on  some  one- 
dimensional  problems  Involving  the  Impinging  of  a  shock 
of  constant  strength  on  an  Interface.   The  method  furnished 
the  transmitted  and  reflected  shocks  fairly  accurately. 

We  contemplate  trying  the  method  for  a  two-dimensional 
flow  involving  shear. 
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Problem  No.  1.19:   On  the  rates  of  convergence  and  the 

"stability  properties  of  finite  difference 
solutions  to  hyperbolic  systems.   (Difference) 

Investigator:      W.  Gersch 

The  objective  of  this  work  was  to  Investigate  the 
characteristics  and  relative  merits  of  the  "centered"  and 
the  "staggered"  finite  difference  schemes  as  they  are  applied 
to  first  order  hyperbolic  systems.   Particular  attention 
was  focused  on  the  rates  of  convergence  and  the  stability 
properties  of  both  schemes. 

The  scheme  referred  to  as  a  centered  scheme  Is  one  In 
which  the  time  and  space  derivatives  are  centered,  the 
staggered  scheme  has  derivatives  centered  In  space  and  staggered 
In  time.   Both  schemes  are  applied  to  rectangular  nets. 

These  schemes  were  applied  systematically  to  the  study 
of  a  one  dimensional  time  dependent  compressible  fluid 
flow  problem.   Machine  calculations  were  made  In  an  attempt 
to  provide  experimental  results  to  supplement  an  analytic 
examination  of  the  finite  difference  schemes. 

A  theoretical  analysis  of  the  application  of  each  of 
the  finite  difference  schemes  to  the  hyperbolic  system 
yields  criteria  for  the  determination  of  the  mesh  ratio  to 
Insure  convergence  and  stability  of  the  numerical  computa- 
tions. In  addition  It  Is  possible  to  extract  a  measure  of 
the  error  of  approximation  of  the  finite  difference  scheme 
to  the  partial  differential  equation. 

The  motivation  for  examining  these  schemes  is  straight- 
forward.  The  convergence  and  stability  of  the  staggered 
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scheme  have  been  analyzed  quite  generally.   The  finite 
difference  approximations  to  the  partial  differential 
equation  for  the  staggered  scheme  are  correct  to  order  of 
the  square  of  the  mesh  width  0(A).   The  centered  scheme 
which  lacks  a  rigorous  general  convergence  proof  is  cor- 
rect to  within  O(A^),  and  therefore  suggests  closer 
approximation  to  the  solution  of  the  partial  differential 
equations.   It  was  thought  desirable  to  get  some  experi- 
mental evidence  on  the  merits  of  each  scheme. 

The  case  examined  experimentally  was  an  Initial 
value  problem  where  the  first  derivatives  do  not  satisfy 
a  Lipschltz  condition.   (This  case  has  as  yet  not  been 
analyzed  successfully) .   An  investigation  of  the  convergence 
and  stability  characteristics  of  finite  difference  schemes 
was  made  and  attention  was  focused  on  the  growth  and 
decay  of  errors  both  within  a  disturbance  region  and  in 
a  region  of  quiet. 

The  experimental  results  suggec-t  that  the  centered 
difference  scheme  does  in  fact  give  a  generally  closer 
approximation  to  the  solution  of  the  partial  differential 
equation  than  the  staggered  scheme.   As  expected,  the 
errors  do  propagate  along  the  characteristics.   Each 
scheme  exhibits  a  particular  behavior  at  the  characteristic 
crossings  and  the  machine  calculations  were  valuable 
here  in  enforcing  a  speculation  concerning  the  random 
Insertion  of  errors. 
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An  additional  Important  observation  was  that  it 
appears  that  there  is  an  optimum  practical  closeness  to 
the  numerical  solution  of  nonlinear  hyperbolic  equations 
whose  derivatives  do  not  satisfy  a  Llpschitz  condition 
and  that  it  is  not  desirable  to  attempt  an  exchange  of 
an  increase  of  computing  time  for  a  closer  numerical 
approximation. 

One  disquieting  result  is  that  neither  scheme 
exhibited  monotonic  convergence  characteristics  within 
the  disturbance  region.   This  suggests  that  additional 
computing  be  done  for  better  behaved  cases  to  provide 
reassurance  that  the  first  impression,  that  the  centered 
scheme  is  in  fact  superior  to  the  staggered  scheme  for 
nonlinear  hyperbolic  systems,  is  on  a  firmer  footing. 
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Problem  No.  1.20:   Numerical  methods  tests 
Investigator:      R.  D.  Rlchtmyer 

Under  this  heading  is  grouped  a  niomber  of  investiga- 
tions which  have  been  going  on  at  various  times  during  the 
past  two  years,  in  part  under  other  numbers  (NYU  1.15 
and  LASL  6. 3).   Most  of  the  results  have  been  presented  at 
various  seminars  and  meetings.   A  few  of  them  will  be 
described  in  a  forthcoming  book  on  Numerical  Methods  for 
Initial-Value  Problems. 

1.20.1  Pseudo-viscosity  method  for  shocks.   When 
finite-difference  methods  are  applied  to  fluid-dynamics 
calculations  using  the  pseudo-viscosity  method  (Von  Neumann 
and  Richtmyer,  Journal  of  Applied  Physics,  1950)  the  ques- 
tions of  accuracy  and  stability  become  Important  and  cannot 
be  completely  settled  by  analytical  methods.   A  code  was 
prepared  for  following  a  one-dimensional  running  wave 
(shock)  in  Lagrangean  coordinates  for  various  values  of 
a)  shock  strength,  b)  time  increment,  c)  ratio  of  specific 
heats,  d)  coefficient  of  pseudo-viscosity.   The  machine 
gives  profiles  of  pressure,  volume,  entropy,  fluid  velocity 
and  obtains  the  shock  speed  and  asymptotic  values  of 
pressure  and  the  like  by  least-squares  and  averaging. 

It  was  found  that  for  reasonable  values  of  the  above 
parameters  the  shock  speed  agrees  with  the  Rankine-Hugoniot 
value  to  a  fraction  of  a  percent;  the  pressure  rise  occurs 
In  two  to  three  space  intervals;  and  behind  the  shock  the 
pressure  oscillates  about  the  correct  value  with  oscilla- 
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tion  amplitudes  of  1  to  5  percent  (depending  on  conditions) 
of  the  pressure  Jump.  The  results  indicate  that  the  method 
should  be  satisfactory  for  many  problems.  « 

Stability  was  investigated  recently,  using  this  code 
to  determine  the  critical  value  of  the  time  increment  for 
various  shock  strengths.   The  results  agree  qualitatively 
with  the  approximate  theoretical  stability  criterion  of  Von 
Neumann  and  Richtmyer,  but  the  equations  are  slightly  less 
stable  for  weak  shocks  and  slightly  more  stable  for  strong 
shocks,  than  was  expected. 

1.20.2   "Imp".   A  family  of  implicit  finite-difference 
equations  was  applied  to  the  solution  of  the  non-linear 
parabolic  equation 

h  da 

where  a  is  a  positive  integer  (=  5  for  the  cases  run  so 
far).   Initial  and  boundary  values  are  taken  from  a  certain 
analytic  solution.   The  machine  prints  the  approximate  and 
exact  solutions  in  parallel  columns  every  ten  cycles.  (The 
exact  solution  is  obtained  by  Newton-Raphson  inversion  of 
a  power  series).   The  agreement  is  very  good,  considering 
the  coarseness  of  the  mesh  used,  and  the  stability  require- 
ment agrees  well  with  the  one  obtained  by  the  usual  heuristic 
arguments. 

1.20.5  Diophantus.   If  i     is  an  irrational  number, 
the  sequence  x  =  fractional  part  of  n^  is  asymptotically 
uniformly  distributed  in  the  internal  0  <  x  <  1.   There  is 
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reason  to  believe  that  for  certain  choices  of  i     the  fluc- 
tuations of  this  distribution  are  smaller  than  would  be  the 
case  for  a  random  sequence.   Because  of  possible  application 
to  model  sampling( "pseudo  Monte  Carlo"),  the  mean-square 
fluctuations  were  computed  for  sequences  x^,  with  n 
running  up  to  several  thousand,  for  various  choices  of  i. 
Double  precision  arithmetic  was  used.   The  fluctuations 
are  noticeably  small  when   ^  is  a  quadratic  or  cubic 
irrational  or  a  rational  linear  function  of  e. 

1.20.4  n^^  root.   A  double-precision  root  extractor 
for  use  with  "diophantus". 

1.20.5  Monte  Carlo  pilot  run.   Early,  preliminary 
version  of  the  resonance  escape  calculation. 

1.20.6  2-dimensional  shock.   (originally  LASL  6.3) 
A  heartbreakingly  unsuccessful  attempt  to  develop  a 

partly  Euleriah,  partly  Lagrangean  method  of  treating  fluid- 
dynamical  problems  in  two  space  variables  and  time. 

1.20.7  Hex  relaxation  convergence  test.   The  over- 
relaxation  or  "extrapolated  Liebman"  method  was  applied  to 
elliptic  problems  of  the  Helmholtz  type  in  two  space 
variables  using  a  hexagonal  net  of  points.   Automatic 
rescallng  as  the  residual  is  quenched  allows  determination 
of  the  fundamental  mode  and  optlmtun  over-relaxation  factor 
with  any  desired  accuracy.   Intended  mainly  for  use  with 
1.20.8  below. 
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1.20.8  2-dlmen3lonal  diffusion.   Various  Implicit 
schemes  of  allegedly  high  accuracy,  using  a  hexagonal  space 
net,  were  applied  to  the  diffusion  equation  in  an  attempt 
to  find  Improved  methods  for  solving  the  multl-dlmenslonal 
age-dlf fusion  equation  In  reactors.   Three  difference  schemed, 
using  21  points,  l4  points,  and  9  points,  respectively,  were 
Investigated  for  a  problem  In  which  the  exact  solution  Is 
known  and  could  be  computed  by  the  machine  (by  power  series) 
for  comparison.   The  Improvement  by  use  of  "high-accuracy" 
formulas  was  disappointingly  small. 

1.20.9  Stability  demonstration.   Mainly  for  Instinac- 
tlon  purposes.   A  simple  heat  flow  problem  Is  solved  by  the 
explicit  equation.   The  time  Increment  can  be  chosen  above 
or  below  the  stability  limit.   The  machine  also  calculates 
the  exact  solution  by  Fourier  series. 

1.20.10  Neutron  transport  integral  equation.  The 
odd  solutions  of  the  eigenvalue  problem 

a 

Tvd)  (s)  =  /  EO|s-r|)(t.(r)dr 
-a 


where 


E(x)  =  r  (e-^)  ^ 

X  ^ 


are  related  tc  the  normal  modes  of  neutron  multiplication 
(or  decay)  in  a  one-group  bare-sphere  problem  with  isotropic 
scattering.   An  Iterative  scheme  was  coded  whereby:  a)  the 
first  two  odd  elgenf unctions  (corresponding  to  the  two 
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highest  values  of  x)  are  found  concurrently,  b)  p  Is 
automatically  adjusted  during  the  calculation  to  make  the 
highest  eigenvalue  converge  to  unity,  c)  In  calculating 
the  fundamental  eigenf unction,  Tchebyscheff  polynomials 
are  used  to  suppress  the  higher  elgenfunctlons .   The 
Integral  Is  evaluated  as  a  Stleltjes  integral  to  avoid 
Inaccuracy  because  of  the  logarithmic  singularity  of  the 
Integrand.   Even  when  a  quite  fine  mesh  is  used,  the  method 
converges  rapidly   The  principal  use  is  to  provide 
accurate  standards  of  comparison  for  the  asymptotic  behavior 
of  the  solutions  of  the  transport  equation  given  by  1.20.11 
below. 

1.20.11   Method  of  characteristics  for  the  time- 
dspendent  transport  equation.   (Continuation  of  work  started 
at  Los  Alamos).   The  transport  equation  for  a  bare  sphere 
of  isotropically  scattering  material  is 

(|  ^  ^  ^  -^  ^^  f  (x,y,t)  =  6f{i^^,t)   , 

where 

r 


y(r,t)  =1^  /  ^ix,/r^~7   ,t)   , 


-r 


2   2     2 
and  is  to  be  solved  in  the  semicircular  region  y  >  0,  x  +y  <  a 

with  the  boundary  conditions 


r(-/a^-y^  >  y>t)  =0       (0  <  y  <  a)  . 

A  direct  integration  method  which  follows  the  characteristics 
(these  are  the  lines  y  =  constant,  x-ct  =  constant)  has 
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been  coded  for  test  of  the  method.   Non-trivial  analytic 
solutions  for  comparison  are  not  knovm,  but  information  can 
be  obtained  by  successive  refinement  of  the  mesh  and  by 
comparison  of  the  asymptotic  form  of  the  solution  with 
functions  obtained  from  1.20.10.   Results  so  far  indicate 
that  the  method  is  accurate  enough  and  fast  enough  that 
it  may  compte  with  other  methods,  such  as  the  S^  method 
of  Carlson,  for  some  purposes. 
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Problem  No.  1.22:   Schrodlnger  wave  equation  (Schemer) 
Investigator:      M.  Salkoff 

The  purpose  of  this  project  Is  to  set  up  a  program 
for  solving  Schrodlnger 's  radial  equation  for  potentials 
that  are  repulsive  near  the  origin.   The  solutions  for 
such  a  potential  are  of  Interest  In  scattering  problems. 
The  problem  under  study  here  Is  that  of  low  energy- 
Inelastic  collisions  between  Hg  molecules.   This  pro- 
blem has  been  examined  by  E.  Bauer  In  research  report 
CX-17,  published  by  the  Institute  of  Mathematical  Sciences. 
In  that  report,  the  WKB  approximation  Is  used  to  solve 
the  radial  equation.   It  Is  felt  that  the  more  accurate 
solutions  obtained  by  the  Unlvac  will  determine  the 
usefulness  of  the  approximation  employed  In  CX-17. 

The  radial  equation  can  be  put  In  the  form 

(1)  ^^  ^   [k^  -  imvirl  .  ij[i±ll]f  (r)  =  0. 

dr*^  h^      r*^ 

2 
Solutions  are  needed  for  various  values  of  k  ,  which 

Is  proportlal  to  the  energy  of  the  incident  Hp  molecule, 
and  I,   which  Is  proportional  to  Its  angular  momentum. 
The  solutions  obtained  are  used  to  calculate  the  Inelastic 
cross-section  by  means  of  simple  perturbation  theory 
and  some  appropriate  perturbation  potential.   This  cross- 
section  will  then  be  compared  with  that  obtained  by  Bauer 
In  CX-17  using  the  WKB  approximation.   The  quantity  of 
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Interest  Is  the  matrix  element  of  the  perturbation  poten- 
tial between  the  initial  and  final  states  of  the  system. 
The  latter  are  in  effect  the  solutions  of  equation  (l) 
for  the  initial  and  final  energies  of  the  incident  particle. 
Since  many  values  of  I     are  needed  in  order  to  calculate 
the  cross-section  at  the  energies  considered  ,  the  Univac 
is  needed  for  the  laborious  series  of  integrations. 

A  Runge-Kutta  integration  scheme  was  used  to  solve 
equation  (l).   The  solutions  so  obtained  were  checked  both 
by  hand  computations  and  by  comparing  with  the  WKB  method 
which  is  fairly  accurate  for  relatively  large  r.   Both 
methods  of  checking  agreed  closely  with  the  Univac  results. 

As  a  further  check,  two  intervals  of  integration  were  used 

2 
to  solve  equation  (l),  for  one  set  of  values  of  k   and 

I.      Comparison  indicated  that  the  solutions  are  good  to 
at  least  four  significant  figures  when  the  solution  is 
of  order  unity. 

The  Univac  program  in  its  present  state  does  not  yet 
give  the  correct  value  of  the  perturbation  matrix  element 
although  the  solutions  of  equation  (l)  do  seem  to  be 
correct,  as  indicated  above.   The  final  total  cross- 
section  calculated  by  means  of  these  matrix  elements  does 
not  vary  correctly  with  energy.   (The  variation  of  cross- 
section  with  energy  is  well  known,  but  the  quantitative 
details  are  not.)   The  error  has  not  yet  been  found. 
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Problem  No.  1.25:   Complex  elgenvalued  differential  equation 

(Friedman)     "" 

Investigators:     S.  Schechter 

This  Is  concerned  with  the  problem  of  finding  a 
numerical  solution  to  the  following  analytic  problem: 

Find  mln  Im  j  >v(w)  \  where  w(t,x)  Is  a  solution 


of 


w 


d  w  ^    (j^2(^)  .  _^  )  „•=  0  (1) 


(») 


dt^  f^ 


lim   |dw  _  i^^i^  J  ^0  (2) 

t->oo  'dt   ^^^o^    ^  ^^ 


\^    w(a,7v)  =  0  (5) 

where  k^(t)  =  Rq  (l+ae'^^^"^0  >  and  a,  k^,  a,  p  are 
given  positive  constants. 

This  problem  arose  out  of  investigations  in  the 
theory  of  tropospheric  refraction  of  radio  waves. 

It  can  be  shown  that  there  exists  a  unique  solution 
to  the  above  problem  without  condition  ()).   We  may 
therefore  regard  the  X.   as  eigenvalues  with  the  pro- 
perty:  w(a,X.)  =  0.   Since  w(t,x)   is  analytic  in  t 
and   A  for  t  >  a  >  0,  these  >,.   are  isolated. 

It  can  also  be  shown  that  if  a  solution  of  (*), 
(w,>,)  exists,  then 

Im 


^^.(e^l'lAoM) 
_  2ako/3 


where  A  =  e 
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One  method  that  has  been  used  in  eigenvalue  problems 
arising  from  two-point  boundary  conditions  is  to  write 
the  equation  in  difference  form  and  solve  a  matrix  pro- 
blem.  As  an  experiment  to  check  this  method  the  problem 
w"  +  A  w  =  0   w(o)  =  w(a)  =  0   was  solved  by  Givens' 
code  for  various  mesh  spacings.   The  results  showed  that 
the  eigenvalues  with  the  smallest  magnitude  yielded  about 
5  digits  accuracy  for  37  mesh  points  but  those  with 
largest  magnitude  had  over  100  percent  error.   Neverthe- 
less one  can  use  several  of  these  values  as  initial  guesses 
for  some  Integration  schemes. 

An  integration  scheme  is  now  being  coded  which  will 
integrate  (l)  from  two  different  starting  points.   One 
point  is  in  a  neighborhood  of  infinity  and  the-  initial 
values  are  obtained  from  an  asymptotic  expression  for  w. 
With  a  choice  for  X  one  would  seek  a  point  t  =  a  such  that 
w(a,>v)  =  0.   The  other  point  of  departure  is  to  start  at 
t  =  a  with  condition  (3)  and  w'(a,>^)  =  c  prescribed.   The 
outward  integration  would  check  for  condition  (2)  in  the 
neighborhood  of  t  =  oo . 
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Problem  No.  1.2^:   Solution  of  integral  equation  (Toff) 
Investigators:     G.  E.  Forsythe,  M.  Salkoff 

In  order  to  solve  the  Integral  equation 

(1)  f(x)  =  g(x)  +  /  K(x,x')  f(x')clx' 

inversion  of  the  integral  operator  by  matrix  methods 
suggests  itself.   This  may  be  a  very  time-consuming  pro- 
cedure, however.   Various  gradient  methods  of  solution  of 
operator  equations  of  the  form 

(2)  V(x)  =  Lu(x) 

have  been  investigated.    The  present  project  was  undertaken 
in  the  hope  that  such  a  method  might  yield  a  solution  of 
equation  (l)  more  rapidly  than  matrix  inversion.   The 
basis  of  such  a  hope  lies  in  the  fact  that  gradient 
methods  of  solution  of  operator  equations  do  not  explicitly 
solve  for  the  inverse  of  the  operator.   Rather,  the  solu- 
tion vector  itself  is  found  by  an  iterative  process,  and 
this  may  entail  fewer  multiplications  and  consequently 
less  time  than  matrix  inversion. 

The  following  iteration  scheme  was  originally  used: 
Suppose 

(3)  V  =  L  u  ;  V  given,  u  sought 

Take 

(4)  V^  =  V 

1.   Forsythe  and  Forsythe,  NBS  report  AMS  39,  Sept.  195^, 
p-  55>  and  bibliography  therein. 
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(5)  V^  =  V^_^   "  ^  ^n      where 

(6)  u^  =  A^  L*^  V^_^       (L*^  =  L  transpose). 
To  determine  X^,  take 

(7)  ^^n'^  ^n^  ~  *^     where  the  brackets  represent 

an  inner  product. 

Using  (5)  for  V^  and  (6)  for  u^,  we  have 


(Vl   -    ^n  ^  L^  Vl'    "^  ^^  Vl) 


0 


or 


,  _  (^Vi'^Vi)       lA 


n-l 


"      (L  l\.,,l  lX.,)      IllViI' 


It  follows   that 
n 


and 


provided 


21  u,  .L-1(V„-  V„) 

1=1 


"  -1  1 

lim  y~  u,  =  L  -^V^  =  L  -^V  =  u 
n->oo  fa  ^       ° 


lim  L"-^  v„  =  0  . 
n.>oo      " 


Two  slight  modifications  were  introduced  in  the  hope  of 
speeding  up  the  convergence  of  the  iteration  procedure: 

(1)  6   -  acceleration  of  the  remainder  vectors  V  ;  and 

(2)  introduction  of  an  arbitrary  multiplicative  factor  p 
in  the  expression  for  x  ,  which  is  varied  from  1  to 
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about  0.8. 

The  results  were  not  very  good,  and  the  method  as  It 
stands  does  not  compare  favorably  with  matrix  inversion. 

Nevertheless,  results  obtained  under  modification  (2) 

2 
above  are  similar  to  those  obtained  by  S.  Stein,   even  to 

the  best  value  of  p. 


2.   Ibid.,  especially  references  (7)  and  (11)  therein, 
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Problem  No.  1.25:   Numerical  methods  for  general  transport 

equations  (Tramp) 

Investigators:     H.  B.  Keller,  B.  Wendroff  (of  L.A.S.L. ) 

Many  problems  concerned  with  weapons,  weapons  effects, 
shielding  and  related  fields  require  the  solution  for  the 
density  0,  of  some  type  of  particles  as  defined  by  the 
transport  equation  in  the  general  form 

(^Mv  7^  "^  ^  *  ?  -^  6(t,r,v)[  <ji{v  ,v  ,v  ,(x>)    =  P(j2f,  t,r,v,(o) . 

Here  r  is  a  position  vector,   v  is  the  speed  of  the 
particles,   m  a  unit  vector  in  the  direction  of  motion 
of  the  particles,  6     is  a  cross-section  for  the  loss  of 
particles,  say  by  fission,  scattering,  emission  or  any 
combination  of  these  processes.   Since  analytic  solutions 
of  such  an  equation  are  known  only  for  extremely  special- 
ized cases,  efficient  numerical  methods  for  the  solution 
of  these  equations  are  of  great  value.   The  present 
program  is  an  attempt  to  investigate,  devise  and  improve 
methods  for  treating  such  problems. 

As  equation  (l)  contains  seven  independent  variables 
numerical  procedures  to  be  used  on  existing  or  proposed 
computers  must  be  specialized  to  less  general  configura- 
tions or  physically  simpler  phenomena  (i.e.  steady  states, 
few  velocities,  weak  angular  dependence,  etc.)   In  partic- 
ular, we  consider  cases  of  spherical  symmetry  or  infinite 
plane  symmetry  in  which  case  the  Independent  variables  are 
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reduced  to  four:  t,  r,  v,  p.,  where  r  Is  radius  or 
distance  from  the  plane  of  origin  and  \i     Is  the  cosine 
of  the  angle  between  position  vector  and  particle  velo- 
city.  A  steady  state  form  of  the  resulting  equations 
(I.e.  0   Independent  of  t)  has  been  studied  theoretically 
and  numerically  In  some  detail  (in  the  case  of  plane 
geometry  where  P  represents  only  scattering  and/or 
emission).   This  program,  called  Achilles,  was  successful 
In  that  a  stable  convergent  numerical  procedure  was 
obtained  and  used  In  experimental  and  practical  problems. 
In  the  present  work  we  consider  spherically  symmetric  one- 
velocity  problems  (i.e.   v  enters  only  as  a  fixed  para- 
meter) with  scattering  or  emission.   The  extension  of 
methods  devised  In  this  study  to  plane  geometries  Is 
Immediate. 

A  basic  approach  to  the  resulting  equation.  In  this 
case 

2 
has  been  presented  by  Bengt  Carlson  .   We  have  modified 

his  formulation  and  are  studying  the  stability  and  con- 


1.  H.  B.  Keller  and  J.  Heller:   On  the  numerical  Inte- 
gration of  the  neutron  transport  equation,  AEC 
Computing  Facility,  Institute  of   Mathematical  Sciences, 
New  York  University,  NYO-6481. 

2.  B.  Carlson:   Solution  of  the  transport  equation  by  S 
approximations,  Los  Alamos  Scientific  Laboratory, 
LA- 1891. 

-  29  - 


vergence  of  the  resulting  scheme.  Numerical  experiments 
are  also  being  conducted  to  determine  the  best  practical 
procedures  and  to  examine  various  features  which  do  not, 
at  present,  seem  amenable  to  theoretical  analysis. 

Our  approach  has  been  to  difference  equation  (2)  in 
the  n-variable  first.   The  resulting  differential-differ- 
ence equations  are  then  seen  to  form  a  hyperbolic  system 
of  first  order  linear  partial  differential  equations,  the 
number  of  dependent  variables  being  equal  to  the  number 
of  ^.-points  taken  in  the  interval  -1  <  p,  <  1.   These 
equations  are  then  differenced,  essentially  as  in  (2) 
but  with  some  Important  simplifications,  such  that  no 
restrictions  are  placed  on  the  mesh  ratio,  Ar/At.   It 
is  shovm  in  the  plane  case  and  in  the  spherical  case 
excluding  the  origin  that  the  difference  scheme  for  the 
hyperbolic  system  is  stable  and  convergent.   Attempts  at 
extending  the  method  and  analysis  to  include  the  origin 
are  in  progress.   Numerical  experiments  show  no  insta- 
bilities at  the  origin  but,  possibly,  some  other  peculiar 
behavior.   The  question  of  convergence  of  solutions  of 
the  difference-differential  equations  to  those  of  the 
transport  equation  (2)  is  more  difficult  and  will  be 
examined  both  theoretically  and  numerically. 

An  important  consequence  of  these  investigations 
has  been  the  development  of  an  unconditionally  stable 
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explicit  difference  scheme  for  hyperbolic  systems  of  first 
order  quasi-linear  partial  differential  equations.   These 
results  will  be  described  in  detail  in  the  near  future. 
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Problem  No.  1.26:   Meteorological  mixed  boundary- Initial 

value  problem  (MeteorJ 

Investigators:     P.  Fife,  E.  Isaacson,  J.  J.  Stoker 


In  a  general  progra^u  for  developing  numerical  methods 
for  solving  hyperbolic  partial  differential  equations 
in  two  or  more  dimensions,  one  of  the  first  difficulties 
that  must  be  overcome  is  the  proper  treatment  of  boundary 
data  in  the  mixed  boundary- initial  value  problem.   We 
have  begun  an  Investigation  of  how  to  treat  the  boundary 
conditions  in  a  meteorological  problem  that  was  formulated 
by  J.  J.  Stoker. 

The  problem  is  to  determine  whether  it  is  possible  to 
predict  the  behavior  of  a  wedge  of  cold  air  on  the  earth's 
surface  by  solving  numerically  an  initial-boundary  value 
problem  involving  three  equations  derived  by  Stoker  as  an 
approximate  description  of  the  situation.   The  three 
equations  are 

u.  +  uu„  +  vu  +  Ah  =  XV 
V      A      y  Jv 

v^  +  uv^  +  vVy  +  Ahy  =  a(-^  u'  -  u) 

h^  +  (uh)^  +  (vh)y  =  0, 

where  u  and  v  are  the  east  and  north  components  of 
the  (cold)  wind  velocity;  h  is  the  height  above  the 
ground  of  the  discontinuity  surface  between  cold  (below) 
and  warm  (above)  air;   x  and  y  are  distances  with 
respect  to  an  orthogonal  coordinate  system  on  a  plane 
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tangent  to  the  earth's  surface,  the  x-axis  running  east, 
the  y-axls  running  north;  p  and  p'   are  densities  of  the 
cold  and  warm  air;  u'   the  warm  air  velocity;  and 

A  =  g(l  -  £^)   , 

X  =  2  o)  sin  (p    , 

where  cd  is  the  angular  velocity  of  the  earth  and  ^    is 
the  latitude.   The  problem  is  somewhat  novel  since  the 
domain  in  which  the  solution  is  sought  is  bounded  on  the 
south  by  a  "free  boundary". 

The  procedure  has  been  to  put  these  equations  in  a 
difference  form  and  carry  out  the  calculations  for  a 
specific  region  of  the  earth's  surface,  bounded  on  the 
east,  west,  and  north  by  fixed  boundaries  and  on  the 
south  by  the  moving  cold  front.   Various  runs  are  being 
made  corresponding  to  several  hours  each  in  actual  time, 
to  determine  the  effect  of  (l)  the  method  used  of 
handling  the  free  boundary;  (2)  the  number  and  type  of 
boundary  conditions  applied  to  the  fixed  boundaries, 
and  (5)  the  initial  data  used. 

So  far  runs  have  been  carried  out  using  two  distinct 
methods  of  treating  the  free  boundary,  the  results  showing 
little  difference.   A  prescription  of  fixed  u,  v,  and  h 
on  the  east  and  west  boundaries  has  been  used  always, 
but  as  for  the  north  boundary,  runs  have  been  made  with 
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(vh)  prescribed  there  as  a  function  of  x  alone,  as  well 
as  runs  prescribing  u,  v,  and  h  as  fixed.   These  both 
have  produced  somewhat  unsatisfactory  results.   A  run  Is 
shortly  to  be  made  In  which  v  =  0  Is  the  only  "artificial" 
boundary  condition  to  be  prescribed  on  the  north  (the 
other  conditions  being  derived  from  the  equations). 

In  the  Initial  conditions  used,  the  deviation  from 
an  equilibrium  pattern  (a  perfect  wedge  with  u  =  const., 
V  =  O)  has  been  localized  In  the  middle  of  the  region. 
This  deviation  has  usually  been  a  ridge  or  a  trough  In  the 
wedge  of  cold  air. 
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Problem  No.  1.27:   Singular  integral  equation  (Ship) 
Investigators:  K.  Grossman,  E.  Isaacson,  A.  Peters,  J.  J.  Stoker 

We  have  begun  a  program  to  develop  methods  that  will 
be  appropriate  for  the  numerical  solution  of  singular 
Integral  equations  of  the  first  kind. 

For  example,  we  want  to  find  the  solution  5>(x,y) 
defined  in  the  domain  D  with  boundary  B,  such  that 

(1)  K(x,y)  =  ^  (^x)jrg(^,n)-H(n-y)5^^(^,-n)dedTi 

[(^x)2+(Ti-y)2] 

where  K(x,y)   is  given  and  the  Integral  is  evaluated  as 
a  Cauchy  principal  value  integral.   That  is, 

^  =  11m   //  ,  where 
D    €— >0  Dg 

Dg  Is  the  subset  of  D  which  is  obtained  by  deleting 

a  circle  of  radius  €  centered  at  the  singular  point  (x,y). 

Owing  to  the  fact  that  equations  such  as  (l)  do  not 
have  unique  solutions,  we  are  able  to  Impose  another 
restriction  on  5'(x,y).   J.  J.  Stoker  and  A.  Peters 
propose  the  boundary  condition  P{B)    =  0  for  their 
physical  problem. 

Prior  to  starting  on  the  two-dimensional  problem,  we 
are  studying  a  one-dimensional  integral  equation,  for 
which  we  are  developing  methods  that  will  be  useful  in 

1.  These  equations  arise  in  the  study  of  the  motion  of 
thin  floating  ships  in  the  work  of  J.  J.  Stoker  and 
A.  Peters.  They  have  suggested  the  problems  we  are 
considering. 
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problems  of  higher  dimension. 

The   equation 

+1 
(2)  K   (x)   =  -  i       i     -^^     ^^     with  the   aiixillary 

relation,    to   ensiire  xxniqueness, 

1 

(2a)      /  f(t)  dt  =  0,  has  been  solved  nximerically. 
-1 

Method  ki     We  tried  the  simplest  finite  difference  approxi- 
mation for  the  integral  in  (2)  based  on  the  representation 
which  exhibits  the  boundary  behavior,   <P(t)  =  -^ — -   •  After 

introducing  t  =  sin  ©,  x  =  sin  y\,      the  integral  equation 
becomes 

(3)      K(.i„  n)  =  -  i    ;  Air.'l^\  ""i^ 


*i 


(3a)      /  ^   (sin  ©)  d©  =  0  . 

With  the  points   Hi  =  "  ^^  "^  ^  T6  *      i  =  1,  2,  . ..,  10, 

lilt        Tt 

and  similarly  ©<  =  "  20"  "^  "^  ife  »  ^®  °^^  express  the 
integral  equation  in  the  form 

10 

lo 


^h)      K  =  -  ^ 


•i}  sin  V|.      , 

jj__ , +  J,   _ i—  +  J, 


sin  ©.  -  sln>^^   ^1  2  cos  r^^^   ^i 

for  i  =1,  2,  ...,  10  and 
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10 
(lj.a)  ^  y^i   =  0»  where  K^  =  K(sin  >|^)  and  \|r.  =  i|r(sin  ©.), 

We  evaluated  the  principal  value  integral  by  expanding  in 
the  interval   ^i  ~  ^  ^  ^  ^  ^1  "•"  f^  »   *°  first  order  in 
(©  -  y\.  )•   In  this  process  the  ten  vinknown  derivatives 
\|(.  =  i|f  (sin  ©.)  are  introduced.  The  eleven  equations  (ij.) 
and  ()+a)  are  supplemented  by  the  nine  compatibility  condi- 
tions 

(Ub)    ^i+i  -  ^i  =2^  t<  *  ^1+1^ 

for  i  =  1,  2,  •••,   9.   Because  of  the  simple  form  of 
equations  (I4.)  and  (i|b)  it  is  easy  to  eliminate  the  unknovm 
derivative  algebraically.  We  then  are  left  with  a  system 
of  ten  linear  equations  for  the  unknowns   t^  ,   i  =  1,  2, 
...,  10. 

A  sample  calculation  was  performed  for  the  function 
K(x)  =  1.   The  resulting  exact  solution  is  given  by 
ip(t)  =  -  ~~ZIIZ       ^^^   th®  accuracy  of  otir  nxomerical 

solution  is  three  significant  figures.   (See  table  below.) 

Method  B:   Since  it  does  not  appear  feasible  to  generalize 

the  finite  difference  method  described  above,  we  are 

examining  the  following  scheme. 

We  let 

(5)     (D(t)  =  --z=:  ^^r.   +  a^t  +  a^t^  +  a„t^]   , 
Vl-f^ 

We  insert  (5)  in  (2)  and  (2a)  and  determine   a  ,  ...»  a^ 
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by  the  least  squares  criterion  that  they  minimize 

+  1  1 

I  (a^,  ...,  a^)  =    /  [mix)   +  i  ;  -|^  dt]2  p(x)dx   , 

-  1  -1 

where  P(x)   is  a  non-negative  weights  Method  B,  if 

successful,  would  be  rather  easily  adaptable  to  problems 

in  two  dimenions. 


Table  of  Finite  Difference  Solution 
for  K(x)  =1 


t 

<?(t)  computed 

<?(t)  exact 

h  = 

sin  (-81) 

+  6.3085133 

+  6.31375131 

*2  = 

sin   (-63) 

+  1.96098235 

+  1. 96261  Oi|9 

*3  = 

sin   (-45) 

+  0.99917162 

+  1.0 

v° 

sin   (-27) 

+  0.50910373 

+  0.50952514.5 

*5  = 

sin   (-9) 

+  0.15825380 

+  o.l5838i|i; 

sin   (9) 

-  0.15825380 

-  0.158381M 

h  = 

sin    (27) 

-  0.50910373 

-  0.509525i|5 

*8  = 

sin   ikS) 

-  0.99917162 

-   1.0 

*9  = 

sin    (63) 

-  1.96098235 

-  1.962610^9 

ho' 

a  sin   (81) 

-  6.3085133 

-  6.31375131 
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Problem  No.  1,28:   Non-linear  buckling  of  spherical  shells 

under  uniform  external  pressure  (Dimples) 

Investigators:   H,  J.  G-reenberg,  H.  B.  Keller,  E.  Reiss, 

M.  Weissner 


It  is  well  knovm  experimentally  that  a  thin  spherical 
shell  becomes  imstable  and  buckles  under  external  lateral 
pressure  in  a  manner  entirely  at  variance  with  the  results 
of  the  linear  buckling  theory  of  elasticity.   Von  Karman 
and  Tsien,  Priedrichs,  Tsien,  and  Altschuler  have  investi- 
gated this  problem  as  one  in  the  non-linear  theory  of 
elasticity  with  small  strains.   In  the  majority  of  this 
work  the  treatments  were  based  on  strain  energy  considera- 
tions.  Although  none  of  these  attempts  is  really  satisfac- 
tory, they  do  give  a  qualitative  picture  of  the  highly 
non-linear  character  of  the  deformation. 

More  recently,  Kaplan  and  Pung  presented  a  theoretical 
and  experimental  treatment  of  this  problem.   They  confined 
their  attention  to  shallow  spherical  shells.   These  shells 
may  be  considered  as  separate  structural  entities,  or  as 
representing  the  region  of  an  entire  sphere  where  the 
buckling  occurs.   The  theoretical  approach  employed  by 
Kaplan  and  Pimg  was  to  solve  the  governing  differential 

equations  a  pair  of  ordinary  second  order  non-linear 

differential  equations  by  expanding  all  of  the  variables 

in  powers  of  the  center  deflection  and  equating  coefficients 
of  equal  power.   This  results  in  a  sequence  of  linear 
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equations.   Unfortunately,  solutions  to  these  equations 
could  only  be  obtained  for  a  small  range  of  parameters. 
However,  the  experimental  results  revealed  an  interest- 
ing behavior  in  the  buckling  loads  and  modes  of  deformation 
with  respect  to  certain  parameters  of  the  shell. 

The  Dimples  Project  is  concerned  with  solving  the 
governing  differential  equations  for  shallow  spherical 
shells  by  power  series  expansions  as  introduced  by 
Priedrichs  and  Stoker  in  their  treatment  of  buckling  of 
circular  plates.  These  equations  depend  upon  two  parameters. 
One  of  these  parameters  is  related  to  the  external  pressure, 
while  the  other  depends  upon  the  dimensions  of  the  shell. 
The  equations  are  solved  for  several  rsmges  of  the  parameters 
under  boundary  conditions  corresponding  to  a  fixed  edge.   It 
is  interesting  to  note  that  when  the  experimental  results  of 
Kaplan  and  Pung  are  interpreted  in  terms  of  these  parameters, 
there  is  an  improved  correlation  of  their  results  as  well  as 
a  reduction  in  the  scatter  of  the  dr.ta. 

The  technique  of  solution  is  to  expand  the  dependent 
variables  in  power  series  in  terms  of  the  independent 
variable  (polar  angle).   Substitution  of  these  series  into 
the  differential  equations  yields  a  recurrence  relation, 
so  that  all  the  coefficients  may  be  expressed  in  terms  of 
the  first  coefficient  in  each  series.   Satisfaction  of  the 
boundary  conditions  at  the  fixed  edge  results  in  non-linear 
algebraic  equations  in  the  coefficients.  Prom  the  recurrence 
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relations  these  equations  depend  only  on  the  first 
coefficients.  These  algebraic  equations  are  then  solved 
by  an  iteration  scheme  based  on  a  Newton-Raphson  technique. 
The  computation  of  coefficients  and  the  iteration  scheme 
are  carried  out  on  the  Univac* 

Thus  far,  the  computations  have  yielded  the  following 
information:   Load  vs,  deflection,  change  in  curvature, 
and  stress,  all  at  the  center  for  the  deformation  up  to 
the  initial  buckling  load.  These  computations  are  for  a 
rather  wide  range  of  the  dimensional  parameter.  The 
results  are  in  excellent  quantitative  (in  many  cases)  as 
well  as  qualitative  agreement  with  the  experimental  results 
of  Kaplan  and  Pung,   In  addition,  other  interesting  and 
tmusual  behaviors  have  been  revealed.   To  obtain  a  complete 
description  of  the  entire  deformation  a  code  is  now  being 
written  to  obtain  deflection,  slope,  and  bending  and 
membrane  stresses  as  functions  of  the  independent  variable 
for  all  values  of  the  parameters.   Furthermore,  an  attempt 
is  under  way  to  obtain  the  post -buckling  behavior  of  the 
shell. 

It  is  hoped  that  the  results  of  these  computations 
will  give  some  indication  regarding  the  location  and 
character  of  whatever  boundary  layer  effects  may  exist. 
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Problem  No.  1,29s  Underwater  explosion  bubble  (Trouble) 
Investigators:     L.  Pulkerson,  H.  B.  Keller 

A  number  of  approximate  theories  of  the  oscillation 
of  an  underwater  explosion  bubble  have  been  presented.   The 
basis  of  most  of  these  are  the  'assumptions  of  incompressi- 
bility  of  the  water,  homogeneity  of  the  gas  bubble,  and 
neglect  of  gravity  and  free  surface  effects.  The  first 
of  the  above  assumptions  has  been  relaxed  in  some  theories 
and  replaced  by  allowing  only  sonic  disturbances  in  the 
water  (as  in  Problem  1.13)»  As  a  verification  of  these 
theories  we  are  now  undertaking  a  direct  numerical  integra- 
tion of  the  original  hydro dynamical  equations  neglecting 
only  gravity  and  the  free  siirface. 

The  exploded  weapon  is  initially  replaced  by  a  homo- 
geneous sphere  of  high  pressure  gas  with  a  polytroplc 
equation  of  state.  The  surrounding  water,  initially 
homogeneous  and  of  infinite  extent,  has  the  usual  poly- 
troplc equation  of  state  and  ambient  pressure  characteristic 
of  any  desired  depth.  The  numerical  integration,  starting 
from  the  above  configuration,  is  intended  to  follow  the 
bubble  radius  through  the  first  few  oscillations.   It  also 
presents  the  successive  shocks  and  rarefactions  sent  into 
the  water  and  bubble. 
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Problem  No.  1,32:   Mapping  problem 
Investigator:      J.  Moser,  P.  Ungar 

This  program  has  as  its  aim  the  investigation  of  an 
area-preserving  mapping  of  a  simple  algebraic  nature  in 
the  large.  We  are  interested  in  studying  the  influence 
of  a  particular  type  of  nonlinearity  which  is  of  importance 
in  problems  connected  with  the  design  of  the  PPAG  synchro- 
tron. 

The  computed  example  is  a  mapping  of  the  simple  form 

o 
X,  =  (x  +  y-^)  cos  a  +  y  sin  a 

y^  =  -  (x  +  y-^)  sin  a  +  y  cos  a   • 

The  problem  is  to  find  regions  of  stability  for  the  mapping. 
This  question  is  solved  in  a  small  neighborhood  of  the 
origin,  but  the  behavior  of  the  image  points  of  the  iterated 
mapping  in  a  large  region  is  completely  unknown. 
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Problem  No.   1.33*      Binary  mixtures    ( Swish) 
Investigator;  E.   L.    Rubin 

One  of   the   fundamental   problems   In  theoretical 
Investigations   of  the   liquid   state   Is   the   determination 
of   the  radial   distribution  function,    the   probability  of 
finding  two  molecules    separated  by  a  distance     r»      The 
problem  has  been  formally  solved    subject   to   certain 
approximations  by  Klrkwood    and   In  a  slightly   different 
but   equivalent  way  by  Born   and  Green.      The  equations   of 
Bom   and   Green  have  been  generalized  to    two   components 
and  this  Unlvac   calculation  Is  being  used  to   investigate 
some    theoretically  predictable  properties   of   the   two 
component   system.      It   is   believed  that   the  roots  of  a 
certain  transcendental   equation  will  yield   the   temperature 
and  density  at  which  condensation  first   takes  place. 
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